//========================================================================================
// WgsSimulator.cc
//----------------------------------------------------------------------------------------
//
// Generate simulated shotgun reads for a whole gnome.
// Models the variations of input size and 
//
//----------------------------------------------------------------------------------------
// $Id: //CBT-depot/CBT++/src/WgsSimulator.cc#2 $
//----------------------------------------------------------------------------------------
//
//========================================================================================
//.......1.........2.........3.........4.........5.........6.........7.........8.........9
//

#include "WgsSimulator.hh"


void
TraceQualityTemplate::GenerateTraceQuality(float errRate,SequenceQuality &Tq,int readSize)
{
  float p20val = 	-10*log10(errRate);
  Tq.resize(readSize);
  for(int i = 0;i<readSize;i++){
    Tq[i] = (int) p20val;
  }
}

//---------------------------------------------------
// * GetNextTraceQuality
//---------------------------------------------------
// 
void
TraceQualityTemplate::GetNextTraceQuality(SequenceQuality &Tq)
{
  int p20count = 0;
  do{

    // If there are more trace qualitites available, grab the next one.
    if ((mTqIdx < mTqVec.size()) && !(mTqVec.size() == 0)){
      Tq = mTqVec[mTqIdx];
      mTqIdx++;
    }else{
      //cout<<"Read new...mTqFileIdx="<<mTqFileIdx<<" TqIdx="<<mTqIdx<<" TqVec.size="<<mTqVec.size()<<endl;
		
      // We've exhausted the ones read in, so now read in another one.  
      FASTAFile FA;
      mTqVec.clear();
      mTqIdx = 0;
		
      // Recycle if we've exhausted the files.   
      if (mTqFileIdx >= mTqFileList.size()) mTqFileIdx = 0;
		
      //cout<<mTqFileList[mTqFileIdx].c_str()<<flush;
		
      izipstream zin;
      string TqFileName = mTqFileList[mTqFileIdx];
      zin.open(TqFileName.c_str());
      //cout<<"Read quality template "<<TqFileName<<endl;		
		
      if (zin.fail()){
	cout<<mTqFileList[mTqFileIdx]<<" file open failed."<<endl;
	exit(1);
      }

      FA.ReadQual(zin,mTqVec);
      zin.close();
      //cout<<" TqVec.size="<<mTqVec.size()<<endl;
      mTqIdx = 0;
      Tq = mTqVec[mTqIdx];
      mTqFileIdx++;
    }
    p20count = Tq.CountPhred20();
  }while(p20count < mMinP20);	 // Ensure that the Tq grabbed meets a minimal quality. 
}


//---------------------------------------------------
// * RandomBase
//---------------------------------------------------
// Returns a random base.
//
//
char 
WgsSimulator::RandomBase()
{
  char rndBase;
  float rnNum = mRNG.Uniform(0,1);  
  if (rnNum < 0.25){
    rndBase = 'A';
  }else{
    if (rnNum < 0.50){
      rndBase = 'C';
    }else{
      if (rnNum < 0.75){
        rndBase = 'G';
      }else{
        rndBase = 'T';
      }
    }
  }
  return(rndBase);
}

//---------------------------------------------------
// * AddReads
//---------------------------------------------------
// 
void 
WgsSimulator::AddReads(string ClonePrefix,
		       int NumClones, 
		       float MeanCloneSize, 
		       float StDevCloneSize,
		       bool bPaired,
		       ostream &fout,ostream &qout)
{
  FASTAFile FA;

  SequenceQuality Tq,FwdReadQuality,RevReadQuality;
  Sequence FwdRead,RevRead,FwdErrRead,RevErrRead;

  for(int i = 0;i < NumClones;i++){

    // Repeat this until we get a "valid" read, where valid means not excessive Ns. 
    int NCount = 0;
    do{

      // Randomly pick a sequence to get the read from, proportional to sequence size. 
      int globalIdx = (int) mRNG.Uniform(0,mTotalSeqSize);
      int bound = 0;
      for(int j = 0;j< mSourceSequences.size();j++){
	bound+= mSourceSequences[j].mBaseSize;
	if (globalIdx < bound){
	  mSeqIdx = j;
	  break;
	}
      }
		
      int CloneStart = (int) mRNG.Uniform(0,mSourceSequences[mSeqIdx].mBaseSize-MeanCloneSize);
      int CloneEnd = CloneStart + (int) mRNG.Gaussian(MeanCloneSize,StDevCloneSize);
		
      if (CloneEnd < CloneStart) CloneEnd = CloneStart+1000;
      if (CloneEnd < 0) CloneEnd = CloneStart+1000;
      if(CloneEnd >= mSourceSequences[mSeqIdx].mBaseSize) CloneEnd = mSourceSequences[mSeqIdx].mBaseSize-1;
				
      // Prescan the region for N's, skipping it if it contains 75 or more N's
      NCount = 0;	 
      for(int k = 0;k< Tq.size();k++){
	int baseIdx = CloneStart+k;
	if (baseIdx < 0) continue;
	if (mSourceSequences[mSeqIdx].Bit2Base(mSourceSequences[mSeqIdx][baseIdx]) == 'N') NCount++;
      }
      // Prescan the reverse region
      for(int l = 0;l < Tq.size();l++){
	int sourceBaseIdx = CloneEnd-Tq.size()-l-1;
	if (sourceBaseIdx < 0) continue;
	if (sourceBaseIdx >= mSourceSequences[mSeqIdx].mBaseSize) continue;
	if (mSourceSequences[mSeqIdx].Bit2Base(mSourceSequences[mSeqIdx][sourceBaseIdx]) == 'N') NCount++;
      }
		
      if (NCount < 75){
		

	if (mErrRate > 0){
	  // Note: This is just a kludge to wedge in a way to generate some 
	  // reads that aren't based on a quality specification.   Here it is 
	  // assumed that we aren't really making clone pairs but just a single read 
	  // and MeanCloneSize and StDevCloneSize are the 
	  int readSize = (int) mRNG.Gaussian(MeanCloneSize,StDevCloneSize);
	  mTqTemplate.GenerateTraceQuality(mErrRate,Tq,readSize);
	}else{
	  // Read trace quality from actual trace quality file.  This determines 
	  // read size and quality. 
	  mTqTemplate.GetNextTraceQuality(Tq);
	}

	mReadNameCount++;
	GenerateFwdReadWithErrors(CloneStart,Tq,FwdRead,FwdReadQuality);
	SetReadName(FwdRead,ClonePrefix,mSourceSequences[mSeqIdx].mName,'U',mReadNameCount);
	SetReadName(FwdReadQuality,ClonePrefix,mSourceSequences[mSeqIdx].mName,'U',mReadNameCount);
	FA.Write(fout,FwdRead);
	FA.WriteQual(qout,FwdReadQuality);
		
	mReadCounter++;
	if (!(mReadCounter % 1000)){
	  float elapsed = EndTiming(mStartTime);
	  cout<<mReadCounter<<" reads. NumClones="<<i<<" ReadNameCount="<<mReadNameCount<<" Time="<<elapsed<<endl;
	}

	if (bPaired){
	  mTqTemplate.GetNextTraceQuality(Tq);
	  GenerateRevReadWithErrors(CloneEnd,Tq,RevRead,RevReadQuality);
	  SetReadName(RevRead,ClonePrefix,mSourceSequences[mSeqIdx].mName,'R',mReadNameCount);
	  SetReadName(RevReadQuality,ClonePrefix,mSourceSequences[mSeqIdx].mName,'R',mReadNameCount);
	  FA.Write(fout,RevRead);
	  FA.WriteQual(qout,RevReadQuality);

	  mReadCounter++;
	  if (!(mReadCounter % 1000)){
	    float elapsed = EndTiming(mStartTime);
	    cout<<mReadCounter<<" reads. NumClones="<<i<<" ReadNameCount="<<mReadNameCount<<" Time="<<elapsed<<endl;
	  }			 
	}
      }else{
	//cout<<"Found region of >= 75 Ns"<<endl;
      }

    }while(NCount>=75);
  }
}

void
WgsSimulator::SetReadName(Sequence &Read,string &prefix,string &suffix,char chem,int count)
{
  Read.mName="";

  char buf[500];
  sprintf(buf,"%s%s%d%c%07d",
          prefix.c_str(),
	  suffix.c_str(),
	  1,
          chem,
          count);
  Read.mName = buf;
}

void
WgsSimulator::SetReadName(SequenceQuality &ReadQual,string &prefix,string &suffix,char chem,int count)
{
  ReadQual.mName="";

  char buf[500];
  sprintf(buf,"%s%s%d%c%07d",
          prefix.c_str(),
	  suffix.c_str(),
	  1,
          chem,
          count);
  ReadQual.mName = buf;
}


void
WgsSimulator::GenerateFwdReadWithErrors(int CloneStart,
					SequenceQuality &Tq,
					Sequence &Read,
					SequenceQuality &ReadQuality)
{
  Read.clear();
  ReadQuality.clear();

  int ErrCnt = 0;

  int tqIdx = 0;
  int baseIdx = CloneStart;
  int insCnt = 0;
  int delCnt = 0;
  for(int i = 0;i< Tq.size();i++){
    if (baseIdx < mSourceSequences[mSeqIdx].mBaseSize){	

      // check for InDel
      if (mRNG.Uniform(0,1) < 0.001){
	// determine if it is an Ins or a Del
	if (mRNG.Uniform(0,1) < 0.5){
	  // Insert
	  insCnt++;
	  Read.push_back(RandomBase());
			 
	  // Knock a little off of the quality of an ins. 
	  int insqual = Tq[tqIdx] - 10;
	  if (insqual < 0) insqual = 0;

	  ReadQuality.push_back(insqual);
	}else{
	  delCnt++;
	  // Delete.  Skip a tq and a base. 
	  baseIdx++;
	}		 
      }else{
	// check for error
	float rate = pow(10,(-(float)Tq[tqIdx]/10.0));
	if (mRNG.Uniform(0,1) < rate){
	  ErrCnt++;
	  Read.push_back(RandomBase());
	  ReadQuality.push_back(Tq[tqIdx]);
	}else{		 
	  // SourceSequence is saved as a Sequence3Bit to minimize space. 
	  // Bit2Base converts the 3-bit pattern to a character. 
	  Read.push_back(mSourceSequences[mSeqIdx].Bit2Base(mSourceSequences[mSeqIdx][baseIdx]));
	  ReadQuality.push_back(Tq[tqIdx]);
	}
      }
    }
	 
    tqIdx++;
    baseIdx++;
  }

  Read.mName = "";
  Read.mDescription = "";

  char buf[500];
  sprintf(buf,"[%d,%d]",CloneStart,CloneStart+Read.size()-1);
  string sbuf = buf;
  Read.mDescription+=sbuf;
  Read.mDescription+=" Fwd";

  Read.mDescription += " Err: ";
  Read.mDescription += itoa(ErrCnt);

  Read.mDescription += " Ins: ";
  Read.mDescription += itoa(insCnt);

  Read.mDescription += " Del: ";
  Read.mDescription += itoa(delCnt);

  Read.mDescription = Read.mDescription + " Size=";
  Read.mDescription = Read.mDescription + itoa(Read.size());

  ReadQuality.mName = Read.mName;
  ReadQuality.mDescription = Read.mDescription;


}

void
WgsSimulator::GenerateRevReadWithErrors(int CloneEnd,
					SequenceQuality &Tq,
					Sequence &Read,
					SequenceQuality &ReadQuality)
{
  Read.clear();
  ReadQuality.clear();

  int ErrCnt = 0;

  // Create the read sequence. 
  for(int i = 0;i< Tq.size();i++){
    int sourceBaseIdx = CloneEnd-Tq.size()+i-1;
    Read.push_back(mSourceSequences[mSeqIdx].Bit2Base(mSourceSequences[mSeqIdx][sourceBaseIdx]));
  }

  // Read is reverse complement relative to fwd read. 
  Read.ReverseComplement();

  int insCnt = 0, delCnt = 0;

  // Now add the errors
  for(int i = 0;i < Tq.size();i++){
    int sourceBaseIdx = CloneEnd-Tq.size()+i-1;
	 
    if ((CloneEnd-sourceBaseIdx) > 0){
      // check for InDel
      if (mRNG.Uniform(0,1) < 0.001){

	// determine if it is an Ins or a Del
	if (mRNG.Uniform(0,1) < 0.5){
	  insCnt++;
			 
	  // Insert
	  Read.push_back(RandomBase());
			 
	  // Knock a little off of the quality of an ins. 
	  int insqual = Tq[i] - 10;
	  if (insqual < 0) insqual = 0;
			 
	  ReadQuality.push_back(insqual);
	}else{
	  delCnt++;
	  // Delete.  Skip a tq and a base. 
	  i++;
	}		 
      }else{
	// check for error
	float rate = pow(10,(-(float)Tq[i]/10.0));
	//		float rate = -10*log10(Tq[i]);  // KJD this has been here and wrong for a long time!
	if (mRNG.Uniform(0,1) < rate){
	  ErrCnt++;
	  Read[i] = RandomBase();
	  ReadQuality.push_back(Tq[i]);
	}else{
	  //Read.push_back(mSourceSequences[mSeqIdx].Bit2Base(mSourceSequences[mSeqIdx](sourceBaseIdx)));
	  ReadQuality.push_back(Tq[i]);
	}
      }
    }
  }
  
  Read.mName = "";
  Read.mDescription = "";
  char buf[500];
  sprintf(buf,"[%d,%d]",CloneEnd-Read.size()+1,CloneEnd);
  string sbuf = buf;
  Read.mDescription+=sbuf;
  Read.mDescription+=" Rev";

  Read.mDescription += " Err:";
  Read.mDescription += itoa(ErrCnt);

  Read.mDescription += " Ins:";
  Read.mDescription += itoa(insCnt);

  Read.mDescription += " Del:";
  Read.mDescription += itoa(delCnt);


  Read.mDescription = Read.mDescription + " Size=";
  Read.mDescription = Read.mDescription + itoa(Read.size());

  ReadQuality.mName = Read.mName;
  ReadQuality.mDescription = Read.mDescription;
}



/*
// The size of the Svec insert
//
//   Uprimer|--------------------|Rprimer
// ========================================
//
int
WgsSimulator::SvecInsertSize() {
return (int)mRNG.Gaussian(mCxtInsertLength,mCxtInsertSdev);
}

// Given a read, add errors and add it to the Simulator
// TODO: write out out to the file, We will not be using it again...
// For now, keep it in memory...
void 
WgsSimulator::AddErrorsAndPushSeq(Sequence &ErrFreeSeq)
{
int readNum;
  
mReads.GrowVec();     // Add another sequence to the read vector
mReadsQual.GrowVec(); // Add another sequence quality to the read vector

readNum=mReads.size()-1;
mReads[readNum].clear();
mReadsQual[readNum].clear();
AddSequencingErrors(ErrFreeSeq,mReads[readNum],mReadsQual[readNum]);

// update counters
mReadNameCnt++;
mTotalPhred20 +=mReadsQual[readNum].CountPhred20();
}

// Make up a Svec insert
void
WgsSimulator::GenInsertEnds(Sequence *sourceSeq, int &start, int &end)
{
int length = 0;
do{
start=(int)mRNG.Uniform(0,sourceSeq->size());
end=start+SvecInsertSize();

// Clamp the end...
if (start >= sourceSeq->size()) {
start = sourceSeq->size()-2;
}
if (end >= sourceSeq->size()) {
end = sourceSeq->size()-1;
} 
length = end-start;
}while(length < 30); // KJD Want to omit pathologically short reads. 
}

int
WgsSimulator::GenReadLength () 
{
int l;
l=(int)mRNG.Gaussian(mCxtReadAvgLength,mCxtReadSdev);
if (mCxtReadMaxLength>0) {
if (l>mCxtReadMaxLength) {
l=mCxtReadMaxLength;
}
}
return l;
}
      
      

// Genadd a random Universal (forward) read
//
//      U----->     <-----U
//    ==============================
//
// These can be on either strand or running either direction.
// The insert is modeled as running in the forward direction.
void
WgsSimulator::GenaddReadUniversal(Sequence *sourceSeq)
{
Sequence read;
int iStart,iEnd;
int rEnd;
char name[200];
string sname;

GenInsertEnds(sourceSeq,iStart,iEnd);
rEnd=iStart+GenReadLength();

if (rEnd > iEnd) {
rEnd=iEnd;
}

sourceSeq->GetSubSequence(read,iStart,rEnd);
SetReadName(read.mDescription,mReadNameCnt,'U');

// The forward reads can be on either strand, going either direction.
if (mRNG.Uniform(0,1)<0.5) {
read.Complement();
read.Reverse();
}

AddErrorsAndPushSeq(read);
}

// Genadd a Universal/reverse read pair.
//
//         U-------->
//    =================================
//                      <-------R
//
// In reality U and R can be flipped end for end.
// But this function makes the U read pointed _always_ to the right.
void
WgsSimulator::GenaddReadUniversalReverse(Sequence *sourceSeq)
{
Sequence read;
int iStart,iEnd;
int freadEnd, rreadStart;
int iNameCnt;
string sname;
int strand;

strand=1;
if (mRNG.Uniform(0,1)<0.5) {
strand=-1;
}

GenInsertEnds(sourceSeq,iStart,iEnd);
iNameCnt=mReadNameCnt;

// Universal
freadEnd=iStart+GenReadLength();
if (freadEnd > iEnd) {
freadEnd=iEnd;
}
sourceSeq->GetSubSequence(read,iStart,freadEnd);
SetReadName(read.mDescription,iNameCnt,'U');

// Can be on either strand...
if (strand==-1) {
//read.Complement(); KJD this is what harley had
read.ReverseComplement(); // isn't it always reverse+Complement????????
}
AddErrorsAndPushSeq(read);

// URError
if (mRNG.Uniform(0,1)<mCxtPartURError) {
// It comes from a different insert...
GenInsertEnds(sourceSeq,iStart,iEnd);
mCountReadURError++;
}

// Reverse
rreadStart=iEnd-GenReadLength();
if (rreadStart<iStart) {
rreadStart=iStart;
}

//
sourceSeq->GetSubSequence(read,rreadStart,iEnd);
SetReadName(read.mDescription,iNameCnt,'R');

// The errors should be running the other way
// KJD This isn't right, is it? It should never be just reverse. 
//read.Reverse();
// This should be the opposite strand of the Universal read
//if (strand==1) {
//  read.Complement();
// }

if (strand == 1){
read.ReverseComplement(); // isn't it always reverse+Complement????????? KJD
}

AddErrorsAndPushSeq(read);
}

//
// These can be on either strand or running either direction.
// The insert is modeled as running in the Universal direction.
void
WgsSimulator::GenaddReadChimeric(Sequence *sourceSeq1,Sequence *sourceSeq2)
{
Sequence read;
Sequence subseq;

int rLength,rSplice;
int size1,s1,e1,size2,s2,e2;
char name[200];
string sname;
string s1desc,s2desc;

rLength=GenReadLength();
rSplice=(int)mRNG.Uniform(0,rLength);

// The read might be shorter than rLength
size1=sourceSeq1->size();
s1=(int)mRNG.Uniform(0,size1);
e1=s1+rSplice;
if (e1>size1) {
e1=size1;
}
size2=sourceSeq2->size();
s2=(int)mRNG.Uniform(0,size2);
e2=s2+(rLength-rSplice);
if (e2>size2) {
e2=size2;
}

read.clear();
sourceSeq1->GetSubSequence(subseq,s1,e1);
read.Concatenate(subseq);
s1desc=subseq.mDescription;
sourceSeq2->GetSubSequence(subseq,s2,e2);
read.Concatenate(subseq);
s2desc=subseq.mDescription;

read.mDescription=s1desc+s2desc;

//
SetReadName(read.mDescription,mReadNameCnt,'U');

// The Universal reads can be on either strand, going either direction.
if (mRNG.Uniform(0,1)<0.5) {
read.Complement();
read.Reverse();
}

AddErrorsAndPushSeq(read);
}

// Set the readname of a read.
void
WgsSimulator::SetReadName(string &readname,int cnt,char chem)
{
char name[500];
char rtemplate[100];
char info[500];
string dir;
string sname;

//
sprintf(name,"%s%d%c%07d",
mReadNamePrefix.c_str(),
mReadNameAttempt,
chem,
cnt);

// Read template
sprintf(rtemplate,"%s%07d",mReadNamePrefix.c_str(),cnt);

// 
if (chem == 'U') {
dir="fwd";
} 
else if (chem == 'R') {
dir="rev";
}
else {
dir="unknown";
}

sprintf(info,"CHEM: %s DYE: big TEMPLATE: %s DIRECTION: %s",
mReadChem.c_str(),
rtemplate,
dir.c_str());

//  if (mReadNameFormat=0) {
//  sname=name+" "+read.mDescription;
//}
//else {
sname=name;
readname=sname+" "+info+" "+readname;
}

// Generate bunches of reads
void 
WgsSimulator::GenaddCntReadsUR(Sequence *sourceSeq,int freadCnt, int rreadCnt)
{
int i;

for (i=0;i<=freadCnt;i++) {
GenaddReadUniversal(sourceSeq);
}
for (i=0;i<=rreadCnt;i++) {
GenaddReadUniversalReverse(sourceSeq);
}
}

void 
WgsSimulator::GenaddCntReads(Sequence *sourceSeq1,
Sequence *sourceSeq2,
int readCnt,
double readSourceRatio  // % from first
)
{
int r;
Sequence *ssrc;
double rkind;

r=0;
while (r<readCnt) {
//
if (mRNG.Uniform(0,1)<readSourceRatio) {
ssrc=sourceSeq1;
}
else {
ssrc=sourceSeq2;
}
// 
rkind=mRNG.Uniform(0,1);
if (rkind<mCxtPartChimeric) {
GenaddReadChimeric(sourceSeq1,sourceSeq2);
mCountReadChimeric++;
r+=1;
}
else if (rkind<mCxtPartU) {
GenaddReadUniversal(ssrc);
mCountReadU++;
r+=1;
}
else { 
GenaddReadUniversalReverse(ssrc);
mCountReadUR++;
r+=2;
}
}
}

void
WgsSimulator::printCxt()
{
printf("Reads are:\n");
printf("   Insert Size = %4d (+/- %3d)\n",mCxtInsertLength,mCxtInsertSdev);
printf("   Read   Size = %4d (+/- %3d)\n",mCxtReadAvgLength,mCxtReadSdev);
if (mCxtReadMaxLength>0) {
// Zero based
printf("          Max  = %4d (if not 0)\n",mCxtReadMaxLength+1); 
}
}

void
WgsSimulator::printCounts()
{
printf("Read counts:\n");
printf("   UR   = %6d \n",mCountReadUR);
printf("   U    = %6d \n",mCountReadU);
printf("Error counts:\n");
printf("   Chimeric   = %6d \n",mCountReadChimeric);
printf("   UR Errors  = %6d \n",mCountReadURError);
}

//------------------------------------------------------------
// Scavenged Code
//------------------------------------------------------------


//---------------------------------------------------
// * AddSequencingErrors
//---------------------------------------------------
// Adds in sequence errors and quality values
// KJD: BUG  This crashes, sometimes, when the sequences are 
// large.   
// 
void 
WgsSimulator::AddSequencingErrors(Sequence &ErrorFreeSeq,
Sequence &SeqOut,
SequenceQuality &QualOut)
{ 
char name[50];
bool MakeDelete = false;

// Remove any existing values (precautionary)
QualOut.clear();
SeqOut.clear();

// Name the sequence
SeqOut.mDescription = ErrorFreeSeq.mDescription;
QualOut.mDescription = ErrorFreeSeq.mDescription; 

int readLen = ErrorFreeSeq.size();
for(int bNum = 0;bNum < readLen;bNum++){
// See if there is an Ins or Del at this location
if(mRNG.Uniform(0,1) < InDelRate(SeqOut.size()-1,readLen)){	 
// Flip a coin to see if it's an Ins or a Del
if (mRNG.Uniform(0,1) >= 0.5){
// Insert
SeqOut.push_back(RandomBase());
int PhredValue = 
CrudePositionDependentPhredValue(SeqOut.size()-1,readLen);
QualOut.push_back(PhredValue);
		  
// Record the insert in the read name	
sprintf(name,".I%d",bNum);
SeqOut.mDescription += name;
mInsCount[SeqOut.size()-1]++; // Record insert distribution
}else{
// Delete
sprintf(name,".D%d",bNum);
SeqOut.mDescription += name;
MakeDelete = true;
mDelCount[SeqOut.size()-1]++; // Record delete distribution
}
}

// For deletes, we skip this base, otherwise add the base to 
// the sequence. 
if (!MakeDelete){
// Add the next base from ErrorFreeSeq
SeqOut.push_back(ErrorFreeSeq[bNum]);

mBaseCount[SeqOut.size()-1]++; // Record base distribution

int PhredValue = 
CrudePositionDependentPhredValue(SeqOut.size()-1,readLen);
QualOut.push_back(PhredValue);
double errRate = pow(10,-((double)PhredValue/10.0));

// Add error if there is one. 
if (mRNG.Uniform(0,1) <= errRate){
SeqOut[SeqOut.size()-1] = RandomBaseError(SeqOut[SeqOut.size()-1]);
mErrCount[SeqOut.size()-1]++; // Record error distribution
}
}
MakeDelete = false;
}
}



//---------------------------------------------------
// * InDelRate
//---------------------------------------------------
// The In/Del rates here are based on the reported errors 
// in:
//
// REF:
//
// Richterich, P. (1998) Estimation of Errors in "Raw" DNA Sequences: 
// A Validation Study.  Genome Research 8:251-259
//
// The frameshift error rate (In/Del combined?) is reported
// to start in the 1% range, falling to < 0.01 % by base 
// 200, and rising, from 300 to 950, up to 6%.  (Dye Terminator. 
// 4%, 0.03%, 6% for dye primer).    
//
// These reported rates are the probability that an ins 
// occured, or a del occured, 
//
// KJD Note: This should be based on chemistry.  
//
// Based on my own (crude?) mathematical modeling, I found that
// for low error rates FrameshiftErrorRate ~ 0.825*InDelRate. 
// This justifies, roughly, equating the two.  
//
float
WgsSimulator::InDelRate(int pos,
int readLen)
{
float fraction = (float)pos/(float)readLen;

float errorRate = 0;
if (fraction < (1.0/38.0)){
errorRate = 0.005; // 0.5%
}else{
if (fraction < (26.0/38.0)){
errorRate = 0.0001; // 0.01%
}else{
if (fraction < (31.0/38.0)){
errorRate = 0.001; // 0.1 %
}else{
if (fraction < (35.0/38.0)){
errorRate = 0.01; // 1%
}else{
errorRate = 0.02; // 2%
}
}
}
}
return(errorRate);
}


//---------------------------------------------------
// * ReadPositionDependentErrorRates
//---------------------------------------------------
//
void
WgsSimulator::ReadPositionDependentErrorRates(string FileName)
{
ifstream fin;
fin.open(FileName.c_str());
if (fin.fail()){
cout<<FileName<<" file open failed"<<endl;
return;
}  
  
double temp;
do{
fin >> temp;
mErrorRate.push_back(temp);
}while(!fin.eof());
}

//---------------------------------------------------
// * PositionDependentErrorRate
//---------------------------------------------------
//
double
WgsSimulator::PositionDependentErrorRate(int pos,
int readLen)
{
// Map the read position for the given read Length into the table
float temp = (float)readLen/(float)mErrorRate.size();
temp = (float)pos/temp;
int idx = (int)(temp -1);

if (idx < mErrorRate.size()){
return(mErrorRate[idx]);	 
}else{
cout<<"Invalid error rate position\n";
return(0);
}
}

//---------------------------------------------------
// * CrudePositionDependentPhredValue
//---------------------------------------------------
// Crudely model the distribution of phred values 
// along the length of a read.  The read is divided 
// into five zones, and phred values in each zone
// are generated according to a gaussian of particular 
// mean and deviation.  
//
//
int 
WgsSimulator::CrudePositionDependentPhredValue(int pos,
int readLen)
{
float fraction = (float)pos/(float)readLen;

float phredValue = 0;
if (fraction < (1.0/38.0)){
phredValue = mRNG.Gaussian(10,5);
}else{
if (fraction < (26.0/38.0)){
phredValue = mRNG.Gaussian(40,10);
}else{
if (fraction < (31.0/38.0)){
phredValue = mRNG.Gaussian(25,10);
}else{
if (fraction < (35.0/38.0)){
phredValue = mRNG.Gaussian(20,10);
}else{
phredValue = mRNG.Gaussian(10,5);
}
}
}
}
if (phredValue < 5) phredValue = 5;
if (phredValue > 51) phredValue = 51;
return((int)phredValue);
}

//---------------------------------------------------
// * RandomBase
//---------------------------------------------------
// Returns a random base.
//
//
char 
WgsSimulator::RandomBase()
{
char rndBase;
float rnNum = mRNG.Uniform(0,1);  
if (rnNum < 0.25){
rndBase = 'A';
}else{
if (rnNum < 0.50){
rndBase = 'C';
}else{
if (rnNum < 0.75){
rndBase = 'G';
}else{
rndBase = 'T';
}
}
}
return(rndBase);
}


//---------------------------------------------------
// * RandomBaseError
//---------------------------------------------------
// Returns a random base (possibly N) that differs 
// from the original base. 
//
//
char 
WgsSimulator::RandomBaseError(char oldBase)
{
char rndBase = oldBase;

float rnNum = mRNG.Uniform(0,1);  
if (rnNum < 0.25){
rndBase = 'A';
}else{
if (rnNum < 0.50){
rndBase = 'C';
}else{
if (rnNum < 0.75){
rndBase = 'G';
}else{
rndBase = 'T';
}
}
}
if (rndBase == oldBase){
rndBase = 'N';
}
return(rndBase);
}

*/

//============================================================
//============================================================
//============================================================




//  //------------------------------------------------------
//  // * PrintNumContigs
//  //------------------------------------------------------
//  // Write the stats for the number of contigs out to a file.   
//  //
//  void
//  JHGStatTable::PrintNumContigs(int longMean,
//                                string FName)
//  {
//    ofstream fout;
//    fout.open(FName.c_str());

//    fout<<"LongReadMean = "<<longMean<<"\n";

//    int y = 0,x = 0;

//    fout<<"\t\t\t";
//    for(x = 0;x < mCols;x++){
//      fout<<(*this)[x][y].mReadSizeMean<<"\t\t\t";
//    }
//    fout<<"\n";
//    for(y = 0;y < mRows;y++){
//      x = 0;
//      float TargetCoverage = 0.5+0.5*y;
//      fout<<TargetCoverage<<"\t";
//      for(x = 0;x<mCols;x++){
//        fout<<(*this)[x][y].mActualCoverage<<"\t";
//        fout<<(*this)[x][y].mPhred20Mean<<"\t";
//        fout<<(*this)[x][y].mContigs<<"\t";
//      }
//      fout<<"\n";
//    }

//    fout<<"-------PhrapTime-----\n";
//    for(y = 0;y < mRows;y++){
//      x = 0;
//      float TargetCoverage = 0.5+0.5*y;
//      fout<<TargetCoverage<<"\t";
//      fout<<(*this)[x][y].mActualCoverage<<"\t";
//      for(x = 0;x<mCols;x++){
//        fout<<(*this)[x][y].mPhrapTime<<"\t\t";
//      }
//      fout<<"\n";
//    }
//  }

//  //------------------------------------------------------
//  // * PrintContigSize
//  //------------------------------------------------------
//  // Write the stats for the size of contigs out to a file.   
//  //
//  void
//  JHGStatTable::PrintContigSize(string FName)
//  {
//    ofstream fout;
//    fout.open(FName.c_str());

//    fout<<"TCover\tACover\tALength\tNumContig\n";

//    int y = 0,x = 0;

//    fout<<"\t\t\t";
//    for(x = 0;x < mCols;x++){
//      fout<<(*this)[x][y].mReadSizeMean<<"\t\t\t";
//    }
//    fout<<"\n";
//    for(y = 0;y < mRows;y++){
//      x = 0;
//      float TargetCoverage = 0.5+0.5*y;
//      fout<<TargetCoverage<<"\t";
//      fout<<(*this)[x][y].mActualCoverage<<"\t";
//      for(x = 0;x<mCols;x++){
//        fout<<(*this)[x][y].mPhred20Mean<<"\t";
//        fout<<(*this)[x][y].mContigSizeMax<<"\t";
//      }
//      fout<<"\n";
//    }

//    fout<<"-------PhrapTime-----\n";
//    for(y = 0;y < mRows;y++){
//      x = 0;
//      float TargetCoverage = 0.5+0.5*y;
//      fout<<TargetCoverage<<"\t";
//      fout<<(*this)[x][y].mActualCoverage<<"\t";
//      for(x = 0;x<mCols;x++){
//        fout<<(*this)[x][y].mPhrapTime<<"\t\t";
//      }
//      fout<<"\n";
//    }
//  }


//  //------------------------------------------------------
//  // * Allocate
//  //------------------------------------------------------
//  //     
//  void
//  JHGStatTable::Allocate(int cols,
//                         int rows)
//  {
//    JHGStatVec tempCol;
//    JHGAsmSimStats tempCell;
  
//    tempCol.reserve(rows);
//    tempCol.insert(tempCol.end(),rows,tempCell);
  
//    this->reserve(cols);
//    this->insert(this->end(),cols,tempCol);

//    mCols = cols;
//    mRows = rows;
//  }


//  //---------------------------------------------------
//  // * SaveErrorStats
//  //---------------------------------------------------
//  //
//  void 
//  WgsSimulator::SaveErrorStats(string FileName)
//  {
//    ofstream fout;
//    fout.open(FileName.c_str());

//    for(int i = 0;i<1000;i++){
//      fout<<i<<"\t"<<mBaseCount[i]<<"\t"<<mErrCount[i]<<"\t"<<mInsCount[i]<<"\t";
//      fout<<mDelCount[i]<<endl;
//    }

//    fout.close();
//  }





//  //---------------------------------------------------
//  // * GetReadEndpoints
//  //---------------------------------------------------
//  //
//  void
//  WgsSimulator::GetReadEndpoints(int ShortMean,
//                                        long &start,
//                                        long &end)
//  {
//    // Start of a read is uniformly distributed over the target molecule
//    start = (int)mRNG.Uniform(0,mSourceSeq.size());
  
//    // The read length is gaussian distributed around mean. Some 
//    // fraction of reads are long reads.  
//    if (mRNG.Uniform(0,1) <= mLongFraction){
//      end = (int)(start + mRNG.Gaussian(mLongMean,mReadStDev));
//    }else{
//      end  = (int) (start + mRNG.Gaussian(ShortMean,mReadStDev));
//    }

//    if (start >= mSourceSeq.size()) start = mSourceSeq.size()-2;
//    if (end >= mSourceSeq.size()) end = mSourceSeq.size()-1; 
//  }


//  //---------------------------------------------------
//  // * GenerateShotgunRead
//  //---------------------------------------------------
//  //
//  void 
//  WgsSimulator::GenerateShotgunRead2(int readNum,
//                                            Sequence &SeqOut,
//                                            int ShortMean,
//                                            long &ReadStart,
//                                            long &ReadEnd)
//  {

//    GetReadEndpoints(ShortMean,ReadStart,ReadEnd);

//    mSourceSeq.GetSubSequence(SeqOut,ReadStart,ReadEnd);

//    char buf[50];
//    sprintf(buf,"R%d.%d",readNum,ReadStart);
//    string temp = buf; // to facilitate concatenation
  
//    string subseqDesc = SeqOut.mDescription;

//    // Name the specified fraction as reverses. 
//    if (mRNG.Uniform(0,1) <= mReverseFraction){
//      SeqOut.mDescription = temp+".r ";
//      SeqOut.mDescription += subseqDesc;
//    }else{
//      SeqOut.mDescription = temp+".s ";
//      SeqOut.mDescription += subseqDesc;
//    }
//  }

//  //---------------------------------------------------
//  // * GenerateShotgunRead
//  //---------------------------------------------------
//  //
//  void 
//  WgsSimulator::GenerateShotgunRead(int readNum,
//                                           Sequence &Seq,
//                                           int ShortMean)
//  {
//    long ReadStart,ReadEnd;

//    GetReadEndpoints(ShortMean,ReadStart,ReadEnd);

//    mSourceSeq.GetSubSequence(Seq,ReadStart,ReadEnd);

//    char buf[50];
//    sprintf(buf,"R%d.%d",readNum,ReadStart);
//    string temp = buf; // to facilitate concatenation
  
//    string subseqDesc = Seq.mDescription;

//    // Name the specified fraction as reverses. 
//    if (mRNG.Uniform(0,1) <= mReverseFraction){
//      Seq.mDescription = temp+".r ";
//      Seq.mDescription += subseqDesc;
//    }else{
//      Seq.mDescription = temp+".s ";
//      Seq.mDescription += subseqDesc;
//    }
//  }

//  //---------------------------------------------------
//  // * AddSimulatedShotgunReads
//  //---------------------------------------------------
//  // Concatenates new shotgun reads to Reads
//  // 
//  void 
//  WgsSimulator::AddSimulatedShotgunReads(float TargetCoverage,
//                                                float TargetPhred20,
//                                                SequenceVector &Reads,
//                                                SequenceQualityVector &ReadQual )
//  {
//    SimulateShotgunReads(TargetCoverage,TargetPhred20,Reads,ReadQual);
//  }

//  //---------------------------------------------------
//  // * SimulateShotgunReads
//  //---------------------------------------------------
//  // 
//  void 
//  WgsSimulator::SimulateShotgunReads(float TargetCoverage,
//                                            float TargetPhred20)
//  {
//    mReads.clear();
//    mReadsQual.clear();
//    SimulateShotgunReads(TargetCoverage,TargetPhred20,mReads,mReadsQual);
//  }



//  //---------------------------------------------------
//  // * SimulateShotgunReads2
//  //---------------------------------------------------
//  // Generates simulated shotgun reads until a specific 
//  // target coverage is reached.  
//  //
//  void 
//  WgsSimulator::SimulateShotgunReads2(float TargetCoverage,
//                                             float TargetPhred20,
//                                             SequenceVector &Reads,
//                                             SequenceQualityVector &ReadsQual,
//                                             vector<int> &StartLocs,
//                                             vector<int> &EndLocs)
//  {
//    // Emperical factor to make the initial ErrorFreeSeq size be such 
//    // that added errors will result in a trace close to TargetPhred20. 
//    int TargetLength = (int)(TargetPhred20/0.788);
 
//    Sequence ErrorFreeSeq;
//    int readNum = Reads.size(); // In case we are concatenating to Reads
//    int readQualNum = ReadsQual.size();
//    if (readNum != readQualNum) cout<<"ERR:"<<readNum<<" "<<readQualNum<<endl;

//    long StartLoc,EndLoc;
//    long TotalPhred20 = 0;
//    float TotalCoverage=0;
//    // Since size of any given read is probabilistic, we simply loop 
//    // until the desired coverage is reached.  
//    while(TotalCoverage<TargetCoverage){
//      Reads.GrowVec();     // Add another sequence to the read vector
//      ReadsQual.GrowVec(); // Add another sequence quality to the read vector

//      // Generating reads is a random process, so there is no guarantee 
//      // that we will get a read that falls in the range of target phred20 
//      // values we specify.  Thus, we loop until we get it.  I'm taking
//      // TargetPhred20 +/- 2% as adequately close.   
//      // KJD BUG This won't work, as written, with mixed read lengths!
//      float Fraction = 1;
//      int phred20 = 0;
//      int trials = 0; // sanity check
//      do{
//        Reads[readNum].clear();
//        ReadsQual[readNum].clear();
//        Reads[readNum].mDescription = "";
//        ErrorFreeSeq.clear();
		
//        GenerateShotgunRead2(readNum,ErrorFreeSeq,TargetLength,
//                             StartLoc,EndLoc);
				
//        //AddSequencingErrors(ErrorFreeSeq,
//        //						  Reads[readNum],
//        //						  ReadsQual[readNum]);
//        phred20 = ErrorFreeSeq.size();
//        Fraction = (float)phred20/(float)TargetPhred20;
//        trials++;
//      }while(((Fraction > 1.5) || (Fraction < 0.8)) && (trials<60));

//      // sanity check monitoring. If this doesn't happen very rarely, 
//      // something is wrong. 
//      if (trials >=60){
//        cout<<"Trials timeout. Fraction = "<<Fraction<<endl;
//      }

//      TotalPhred20 += phred20;
//      TotalCoverage = (float)TotalPhred20/(float)mSourceSeq.size();
//      readNum++;
//      StartLocs.push_back(StartLoc);
//      EndLocs.push_back(EndLoc);
//    }
//  }


//  void 
//  WgsSimulator::AddSequencingErrorsWithBounds(Sequence &ErrorFreeSeq,
//                                                     Sequence &SeqOut,
//                                                     SequenceQuality &QualOut,
//                                                     float TargetPhred20,
//                                                     float bounds
//                                                     )
//  {
//    int phred20;
//    float Fraction;
//    int trials=0;

//    do {
//      AddSequencingErrors(ErrorFreeSeq,
//                          SeqOut,
//                          QualOut
//                          );
    
//      phred20 = QualOut.CountPhred20();
//      Fraction = (float)phred20/(float)TargetPhred20;
//      trials++;
//    } while ((fabs(1-Fraction)>bounds) && (trials<20));

//    // sanity check monitoring. If this doesn't happen very rarely, 
//    // something is wrong. 
//    if (trials >=20){
//      cout<<"Trials timeout. Fraction = "<<Fraction<<endl;
//    }
//    //cout<<"trials ="<<trials<<endl; // forgot =0
  
//  }


//  //---------------------------------------------------
//  // * SimulateShotgunReads
//  //---------------------------------------------------
//  // Generates simulated shotgun reads until a specific 
//  // target coverage is reached.  
//  //
//  void 
//  WgsSimulator::SimulateShotgunReads(float TargetCoverage,
//                                            float TargetPhred20,
//                                            SequenceVector &Reads,
//                                            SequenceQualityVector &ReadsQual)
//  {
//    // Emperical factor to make the initial ErrorFreeSeq size be such 
//    // that added errors will result in a trace close to TargetPhred20. 
//    int TargetLength = (int)(TargetPhred20/0.788);
  
//    Sequence ErrorFreeSeq;
//    int readNum = Reads.size(); // In case we are concatenating to Reads
//    int readQualNum = ReadsQual.size();
//    if (readNum != readQualNum) cout<<"ERR:"<<readNum<<" "<<readQualNum<<endl;


//    long TotalPhred20 = 0;
//    float TotalCoverage=0;
//    // Since size of any given read is probabilistic, we simply loop 
//    // until the desired coverage is reached.  
//    while(TotalCoverage<TargetCoverage){
//      Reads.GrowVec();     // Add another sequence to the read vector
//      ReadsQual.GrowVec(); // Add another sequence quality to the read vector

//      // Generating reads is a random process, so there is no guarantee 
//      // that we will get a read that falls in the range of target phred20 
//      // values we specify.  Thus, we loop until we get it.  I'm taking
//      // TargetPhred20 +/- 2% as adequately close.   
//      // KJD BUG This won't work, as written, with mixed read lengths!
//      float Fraction = 1;
//      int phred20 = 0;

//      Reads[readNum].clear();
//      ReadsQual[readNum].clear();
//      Reads[readNum].mDescription = "";
//      ErrorFreeSeq.clear();

//      GenerateShotgunRead(readNum,ErrorFreeSeq,TargetLength);

//      AddSequencingErrorsWithBounds(ErrorFreeSeq,
//                                    Reads[readNum],
//                                    ReadsQual[readNum],
//                                    TargetPhred20,
//                                    0.1
//                                    );

//      TotalPhred20 += ReadsQual[readNum].CountPhred20();
//      TotalCoverage = (float)TotalPhred20/(float)mSourceSeq.size();
//      readNum++;			
//    }
//  }




//  //---------------------------------------------------
//  // * SimulateShotgunReads
//  //---------------------------------------------------
//  // Generates a specific number of simulated shotgun reads.  
//  //
//  void
//  WgsSimulator::SimulateShotgunReads(int numReads,
//                                            float TargetPhred20)
//  {
//    SimulateShotgunReads(numReads,TargetPhred20,mReads,mReadsQual);
//  }


//  //---------------------------------------------------
//  // * SimulateShotgunReads
//  //---------------------------------------------------
//  // Generates a specific number of simulated shotgun reads.  
//  //
//  void
//  WgsSimulator::SimulateShotgunReads(int numReads,
//                                            float TargetPhred20,
//                                            SequenceVector &Reads,
//                                            SequenceQualityVector &ReadsQual)
//  {
//    int TargetLength = (int)(TargetPhred20/0.788);
//    Sequence ErrorFreeSeq;
//    for(int readNum = 0;readNum<numReads;readNum++){
//      Reads.GrowVec();
//      ReadsQual.GrowVec();
//      Reads[readNum].clear();
//      ReadsQual[readNum].clear();
//      Reads[readNum].mDescription = "";
//      ErrorFreeSeq.clear();
//      GenerateShotgunRead(readNum,ErrorFreeSeq,TargetLength);
//      AddSequencingErrors(ErrorFreeSeq,Reads[readNum],ReadsQual[readNum]);
//    }
//  }

//  //---------------------------------------------------
//  // * SimulateErrorFreeShotgunReads
//  //---------------------------------------------------
//  // Generates a specific number of simulated shotgun reads.  
//  //
//  void
//  WgsSimulator::SimulateErrorFreeShotgunReads(int numReads,
//                                                     int TargetLength)
//  {
//    mRNG.SeedWithClock();

//    Sequence ErrorFreeSeq;
//    for(int readNum = 0;readNum<numReads;readNum++){
//      mReads.GrowVec();
//      mReadsQual.GrowVec();
//      mReads[readNum].clear();
//      mReadsQual[readNum].clear();
//      mReads[readNum].mDescription = "";
//      ErrorFreeSeq.clear();
//      GenerateShotgunRead(readNum,ErrorFreeSeq,TargetLength);
//      mReads[readNum] = ErrorFreeSeq;
//    }
//  }

//  //---------------------------------------------------
//  // * RunSimulation
//  //---------------------------------------------------
//  //
//  void
//  WgsSimulator::RunSimulation()
//  {
//    // Seed random number generator. 
//    mRNG.SeedWithClock();

//    FASTAFile InFile,OutFile;
//    SequenceVector InSeqVec;

//    int NumCoverageConditions = 16;
//    int NumLengthConditions = 8;

//    mStats.Allocate(NumLengthConditions,NumCoverageConditions);
  
//    // For various read lengths
//    for(int LengthIdx = 0;LengthIdx < NumLengthConditions;LengthIdx++){
//      float TargetPhred20 = 300 + LengthIdx*100;
	 
//      // And various amounts of coverage
//      for(int CoverageIdx = 0;CoverageIdx < NumCoverageConditions;CoverageIdx++){

//        // Delete at beginning so as to leave last constructed set on disk
//        // for examination/debug. 
//        system("rm -f phrap*");

//        float TargetCoverage = 0.5 + 0.5*CoverageIdx;

//        cout<<"Simulating shotgun reads ("<<TargetCoverage<<","<<TargetPhred20;
//        cout<<")..."<<flush;
//        SimulateShotgunReads(TargetCoverage,TargetPhred20);
//        cout<<"done."<<endl;

//        // Record information about the simulated project...
//        mStats[LengthIdx][CoverageIdx].mReadSizeMean = (int)TargetPhred20;
//        //mStats[LengthIdx][CoverageIdx].mNumReads = NumReads;
//        long phred20Total = mReadsQual.Phred20Total();
//        float ActualCoverage = phred20Total/(float)mSourceSeq.size();
//        mStats[LengthIdx][CoverageIdx].mActualCoverage = ActualCoverage;
//        float phred20Mean,phred20Variance;
//        mReadsQual.Phred20Stats(phred20Mean,phred20Variance);
//        mStats[LengthIdx][CoverageIdx].mPhred20Mean = phred20Mean;

//        // Save some stats about the number and size of errors. 
//        string ErrStatName = "err.stat";
//        char temp[50];
//        sprintf(temp,"%f",phred20Mean);
//        ErrStatName += temp;
//        ErrStatName += ".";
//        sprintf(temp,"%f",ActualCoverage);
//        ErrStatName += temp;
//        SaveErrorStats(ErrStatName);

//        // Save simulated reads...
//        cout<<"Writing seq and qual files.."<<flush;
//        OutFile.Write("phrap.seq",mReads);
//        OutFile.WriteQual("phrap.seq.qual",mReadsQual);
//        cout<<"done."<<endl;

//        // Run Phrap on simulated shotgun project...
//        cout<<"=============== Run Phrap  =============="<<endl<<endl;
//        long StartTime;
//        StartTiming(StartTime);
//        system(mPhrapCommand.c_str()); // run phrap on sequences. 
//        float PhrapTime = EndTiming(StartTime);
//        cout<<"=============== Phrap Done =============="<<endl<<endl;

//        cout<<"Reading assembly results.."<<flush;
//        // If phrap fails for some reason, there won't be a contigs file, 
//        // so we just skip and go on. 
//        if (InFile.Read("phrap.seq.contigs",InSeqVec)){
//          cout<<"done"<<endl;
		
//          // Record information about resulting assembly...
//          mStats[LengthIdx][CoverageIdx].mContigs = InSeqVec.size();

//          // Record size of largest contig, which is the last in the 
//          // phrap.seq.contigs file. 
//          int sizeMax = InSeqVec[InSeqVec.size()-1].size();
//          mStats[LengthIdx][CoverageIdx].mContigSizeMax = sizeMax;

//        }else{
//          // Record the failure.
//          mStats[LengthIdx][CoverageIdx].mContigs = -1;
//        }
//        mStats[LengthIdx][CoverageIdx].mPhrapTime = (int) PhrapTime;

//        // Save a table of numContigs vs. num reads and read length
//        mStats.PrintNumContigs(mLongMean,"stats.numcontigs");
//        // Save a table of contig size  vs. num reads and read length
//        mStats.PrintContigSize("stats.contigsize");
//      }
//    }
//  }


// $Log: WgsSimulator.cc,v $
// Revision 1.11  2000/05/01 14:59:12  harley
// Added max read length option.
//
// Revision 1.10  2000/04/13 15:17:16  havlak
// Moved "MakeDelete" flag reset outside of loop test, to fix problem
// where reads were truncated at first deletion error.
//
// Revision 1.9  2000/04/10 23:02:52  harley
// Have to ask donna about this.
//
// Revision 1.8  2000/04/10 15:13:35  harley
// *** empty log message ***
//
// Revision 1.7  2000/04/05 19:44:42  harley
// Changed readnames.
//
// Revision 1.6  2000/04/05 19:08:37  harley
// *** empty log message ***
//
// Revision 1.5  2000/03/31 23:16:16  harley
// VersOne
//
// Revision 1.4  2000/03/31 21:09:18  harley
// *** empty log message ***
//
// Revision 1.3  2000/03/31 19:53:54  harley
// *** empty log message ***
//
// Revision 1.2  2000/03/31 00:38:14  harley
// *** empty log message ***
//
// Revision 1.1  2000/03/30 15:50:43  harley
// Initial revision
//


